Pohlsches Drehpendel mit Zusatzmasse

Differentialgleichung
---------------------

(J0+m*R^2)*phi''+kdG*phi+kdI*I^2*phi'+D*phi-m*g*R*sin(phi)=M0*sin(omega0*t)

->(I) phi'  = omega
  (II)omega'= -kdG*phi-kdI*I^2*phi'-D*phi+m*g*R*sin(phi)+M0*sin(omega0*t)/(J0+m*R^2)

Anfangsbedinungen
-----------------
Die Anfangsbedingungen sind hier auf phi(t=0)=0.5 und phi'(t=0)=0 festgelegt,
da nur das eingeschwungene Verhalten betrachtet werden soll.

Parameter
---------

J0: Trägheitsmoment des Drehpendels
m : Zusatzmasse
R : Entfernung der Zusatzmasse zur Drehachse
D : Federkonstante
kdG: Reibungskoeffizient der Gleitreibung
kdI: Proportionalitätskonstante zwischen I^2 und der von der Wirbelstrombremse verursachten Dämpfung
I: Dämpfungsstromstärke
M0: Amplitude des antreibenden Drehmoments
omega0: Kreisfrequenz des antreibenden Drehmoments 

pohlsches_rad(J0,m,R,D,kdG,kdI,M0,T0,I):=block([dt,lsg,Jg,omega0],

    load(draw),     /* Laden des Paketes 'draw'*/
    load(dynamics), /* Laden des Paktes 'dynamics', welches das Runge-Kutta-Verfahren 4-ter Ordnung beinhaltet.*/
    numer:true,     /* Berechnungen auf numerisch stellen */

    dt:0.03,        /* Festlegen des Zeitschritts*/

    g:9.81,
    omega0:2*%pi/T0,
    Jg:J0+m*R^2,
    c1:D/Jg,
    c2:m*g*R/Jg,
    c3a:kdI*I^2/Jg,
    c3b:kdG/Jg,
    c4:M0/Jg,

    /*Runge-Kutta-Verfahren 4-ter Ordnung, wobei 70 Perioden betrachtet werden */
    lsg:rk([-c1*phi+c2*sin(phi)-c3a*omega-c3b/abs(omega)*omega+c4*sin(omega0*t),omega],[omega,phi],[0.00000001,0.5],[t,0,70*T0,dt]), /*round(70*T0/dt)*/
    lsg:rest(lsg,round(50*T0/dt)),/*round(50*T0/dt)*/

    /* Funktion, welche die Lösung von der Form [t,omega,phi] in die Form [t,phi] überführt */
    f(l):=block([m],
        m:[l[1],l[3]],
        return(m)),

    /* Funktion, welche die Lösung von der Form [t,omega,phi] in die Form [t,omega] überführt */
    w(l):=block([m],
        m:[l[1],l[2]],
        return(m)),

    /* Funktion, welche die Lösung von der Form [t,omega,phi] in die Form [omega,phi] überführt */
    u(l):=block([m],
  m:[l[3],l[2]],
        return(m)),

    /* Funktion, welche die Lösung von der Form [t,omega,phi] in die Form [t,MA] überführt */
    h(l):=block([m],
        m:[l[1],M0*sin(omega0*l[1])],
        return(m)),

    /* Anwendung der Funktionen f, g, h ,w auf lsg */  
    zeit_phi:map(f,lsg),
    zeit_omega:map(w,lsg),
    omega_phi:map(u,lsg),
    zeit_MA:map(h,lsg),

    /* Grafische Ausgabe über das Paket 'draw' */
    wxdraw2d(
        title="phi(t)-Diagramm",
        point_type=dot, 
        xlabel="t/s",
ylabel="phi/rad",
        color=blue,
        points_joined=true,
        points(zeit_phi),
        color=red),
   wxdraw2d(
        title="Phasendiagramm",
        xlabel="phi/rad",
        ylabel="omega/ rad s^-1",
        point_type=dot,
        points_joined=true,
        color=blue,
        points(omega_phi)),
    wxdraw2d(
        title="omega(t)-Diagramm",
        point_type=dot, 
        xlabel="t/s",
        ylabel="omega/rad s^-1",
        color=red, 
        points_joined=true,
        points(zeit_omega),
        color=red))$

      pohlsches_rad(J0,m,R,D,kdG,kdI,M0,T0,I)

pohlsches_rad(0.0021,0.03,0.090,0.0254,0.000063,0.0052,0.002,3.9,0.45);